RepSIRdi <-function(n=1000,d=20,npopV=c(100,100,100),R0=3,I0=c(0,0,0),pmfdi=c(0.2,0.5,0.3),pmfdexc=c(0.2,0.5,0.3),excV=c(1,2,3),facpasym=1/2,propsympt=2/3,pmf1day=c(0.2,0.3,0.5),pmfsend=c(1,1,1,0.6,0.7,1,1,1),pmfstayhome=c(0.1,0.2,0.1),pimm=c(0.2,0.4,0.3)){
  nc<-length(npopV)               
  dimax<-length(pmfdi)
  dvmax<-length(pmfdexc)    
  pinfV<-numeric(nc)              
  excm<-matrix(0,nc,nc)     
  excm[lower.tri(excm,diag=FALSE)]<-excV 
  excm<-t(excm)
  excm[lower.tri(excm)]<-t(excm)[lower.tri(excm)] 
  ratematrix<-function(excm,npopV,nc){
  #This function returns the exchange rate matrix between cities
    rate<-matrix(0,nc,nc)
    for(r in 1:nc){
     rate[r,]<-excm[r,]/npopV[r]
      }
return(rate)
}
pexc<-ratematrix(excm,npopV,nc)
expectdi <- sum((1:dimax)*pmfdi) 
pinfV<-R0/(expectdi*npopV*(propsympt+(1-propsympt)*facpasym)) 
if(length(I0)==1){
 I0V<-rep(0,nc)
 I0V[sample(1:nc,1)]<-I0 }
else{
 I0V<-I0}
 resm<-array(NA,dim=c(n,d,length(npopV),2)) 
 resmIR<-array(NA,dim=c(n,d,length(npopV),2))  
OneSim<-function(){
sa<-array(0,dim=c((2+dimax*2),1+nc*dvmax,nc))   
 imm<-numeric(nc)
 for(p in 1:nc) {
 sa[1:3,1,p] <- c(npopV[p]-I0V[p],0,I0V[p])     
 sa[4:(dimax*2+2),1,p]<-0
 sa[1:(dimax*2+2),(p-1)*dvmax+(1:dvmax)+1,p]<-NA}
  for(f in 1:nc){
  imm[f]<-rbinom( 1, size=sa[1,1,f], prob=pimm[f])
  sa[2,1,f]<-sa[2,1,f]+imm[f]                
  sa[1,1,f]<-sa[1,1,f]-imm[f]}
   resI <- array(0,dim=c(nc,d,2))    
   resIR<- array(0,dim=c(nc,d,2))
for(i in 1:d)
nw<-generatenewlyinfected(sa,dvmax=dvmax,nc=nc,dimax=dimax,pmfdi=pmfdi,pinfV=pinfV,facpasym=facpasym)
nwsymp<-array(0,dim=c(1,1+(nc-1)*dvmax,nc)) 
for(y in 1:nc){
nwsymp[1,,y]<-rbinom(1+(nc-1)*dvmax,size=nw[1,,y],prob=propsympt) } 
nwasymp<-nw-nwsymp             
sa[2,,]<-sa[2,,]+sa[(dimax+2),,]+sa[(dimax*2+2),,]       
sa[5:(dimax+2),,]<-sa[4:(dimax+2-1),,]
sa[(5+dimax):(dimax*2+2),,]<-sa[(4+dimax):(dimax*2+1),,]
sa[4,,]<-0
sa[4+dimax,,]<-0
 after1dayS<-array(0,dim=c(dimax,1+(nc-1)*dvmax,nc))
after1dayA<-array(0,dim=c(dimax,1+(nc-1)*dvmax,nc))
 for(h in 1:nc){
 after1dayS[,,h]<-rmultinomv(n=1+(nc-1)*dvmax,size=sa[3,-((h-1)*dvmax+(1:dvmax)+1),h],prob=pmf1day)
 after1dayA[,,h]<-rmultinomv(n=1+(nc-1)*dvmax,size=sa[3+dimax,-((h-1)*dvmax+(1:dvmax)+1),h],prob=pmf1day) }
for(k in 1:nc){
 sa[2,-((k-1)*dvmax+(1:dvmax)+1),k]<-sa[2,-((k-  1)*dvmax+(1:dvmax)+1),k]+after1dayS[1,,k]+after1dayA[1,,k]
sa[4:(dimax+2),-((k-1)*dvmax+(1:dvmax)+1),k]<-sa[4:(dimax+2),-((k-1)*dvmax+(1:dvmax)+1),k]+after1dayS[(2:dimax),,k]
 sa[(4+dimax):(dimax*2+2),-((k-1)*dvmax+(1:dvmax)+1),k]<-sa[(4+dimax):(dimax*2+2),-((k-1)*dvmax+(1:dvmax)+1),k]+after1dayA[(2:dimax),,k]
 sa[3,,k]<-0
sa[3+dimax,,k]<-0}
for(g in 1:nc){
sa[3,(g-1)*dvmax+(1:dvmax)+1,g]<-NA
sa[3+dimax,(g-1)*dvmax+(1:dvmax)+1,g]<-NA}                            
 if(i<=50){extrainfected<-rbinom(1,size=100,prob=0.01
 else {extrainfected<-0}
 for(h in 1:nc){
 sa[3,-((h-1)*dvmax+(1:dvmax)+1),h]<-sa[3,-((h-1)*dvmax+(1:dvmax)+1),h]+nwsymp[1,,h]  
sa[(3+dimax),-((h-1)*dvmax+(1:dvmax)+1),h]<-sa[(3+dimax),-((h-1)*dvmax+(1:dvmax)+1),h]+nwasymp[1,,h]        
 sa[1,-((h-1)*dvmax+(1:dvmax)+1),h]<- sa[1,-((h-1)*dvmax+(1:dvmax)+1),h]-nwsymp[1,,h]-nwasymp[1,,h]   }        
 sg<-array(0,dim=c(2+dimax*2,nc*dvmax+1,nc))
 for(s in 1:nc){
 sg[1:(dimax*2+2),(s-1)*dvmax+(1:dvmax)+1,s]<-NA }
for(f in 1:nc){
 for(z in 1:(2+dimax*2)){
 sg[z,-((f-1)*dvmax+(1:dvmax)+1),f]<-rbinom((nc-1)*dvmax+1,size=sa[z,-((f-1)*dvmax+(1:dvmax)+1),f],prob=pmfsend[z]) }}
sh<-array(0,dim=c(dimax,nc*dvmax+1,nc))
 for(j in 1:nc){
 sh[1:dimax,(j-1)*dvmax+(1:dvmax)+1,j]<-NA} 
 for(p in 1:nc){
 for(k in 1:dimax){ 
 sh[k,-((p-1)*dvmax+(1:dvmax)+1),p]<-rbinom((nc-1)*dvmax+1,size=sg[k+2,-((p-1)*dvmax+(1:dvmax)+1),p],prob=pmfstayhome[p])}}
sg[3:(dimax+2),,]<-  sg[3:(dimax+2),,]-sh[1:dimax,,]
sa<-send(dvmax=dvmax,sg,pexc=pexc,pmfdexc=pmfdexc,nc)     
 rst<-array(0,dim=c(nc,d,2))
 rstR<-array(0,dim=c(nc,d,2))
 for(rs in 1:nc){
  rst[rs,,1]<-sum(sa[3:(dimax+2),,rs],na.rm=TRUE)             
  rst[rs,,2]<-sum(sa[3:(2*dimax+2),,rs],na.rm=TRUE)                
  rstR[rs,,1]<-sum(sa[2:(dimax+2),1,rs])+sum(sa[2:(dimax+2),(rs-1)*dvmax+(1:dvmax)+1,],na.rm=TRUE)
  rstR[rs,,2]<-sum(sa[2:(2*dimax+2),1,rs])+sum(sa[2:(2*dimax+2),(rs-1)*dvmax+(1:dvmax)+1,],na.rm=TRUE)}       
resI[,i,1]<-rst[,i,1]
 resI[,i,2]<-rst[,i,2]
 resIR[,i,1]<-rstR[,i,1]
 resIR[,i,2]<-rstR[,i,2]}
 return(list(I=resI,IR=resIR,IM=imm)) }
 for(i in 1:n){
 inf<-OneSim()
 mm<-inf$IM
 for(y in 1:2){
 for(h in 1:length(npopV)){
 resm[i,,h,y]<- inf$I[h,1:d,y]
 resmIR[i,,h,y]<- inf$IR[h,1:d,y]}}}
 return(list(I=resm,IR=resmIR,IMM=mm))
generatenewlyinfected<-function(sa,dvmax=3,nc=3,dimax=3,pmfdi=c(0.2,0.3,0.5),pinfV=c(0.001,0.002,0.003),facpasym=1/2){ 
 #This function calculates and returns the newly infected cases for each city 
newinf<-array(0,dim=c(1,1+(nc-1)*dvmax,nc)) 
for(x in 1:nc){
pinfection <- 1-(1-pinfV[x])^sum(sa[3:(dimax*2+2),,x],na.rm=TRUE)*(1- pinfV[x]*facpasym)^sum(sa[(3+dimax):(2+2*dimax),,x],na.rm=TRUE)
newinf[,,x] <- rbinom( (1+(nc-1)*dvmax), size=sa[1,-((x-1)*dvmax+(1:dvmax)+1),x], prob=pinfection)
}  return(newinf)   }
probforsend<-function(nc,pmfdexc,pexc){  
  #By the help of this function, probability values of staying in other cities for a length of time are calculated by using pmfdexc and pexc inputs as described in the general function of the code 
 dvmax<-length(pmfdexc)
 probarray<-array(0,dim=c(dvmax,nc+1,nc))  
 for(p in 1:nc){
 probarray[1:(dvmax),p,p]<-NA }
 for(i in 1:(nc-1)){
 for(j in (i+1):nc){
 probarray[,j,i]<-pmfdexc*pexc[i,j]
probarray[,i,j]<-pmfdexc*pexc[j,i] } }
 for(k in 1:nc) { 
 probarray[1,nc+1,k]<-1-sum(probarray[,-k,k])
 probarray[2:dvmax,nc+1,k]<-NA } 
 return(probarray)}
send<-function(dvmax,sa,nc,pmfdexc,pexc){ 
  #This function returns the new sa array after sending people among the cities
 prb<-probforsend(nc=nc,pmfdexc=pmfdexc,pexc=pexc)       
 trv<-array(0,dim=c(length(sa[,(nc*dvmax+1),1]),nc*dvmax+1,nc))  
for(p in 1:nc){
trv[,(p-1)*dvmax+(1:dvmax),p]<-NA }  
for(t in 1:nc){
 tr<-generatesend(sa[,1,t],prb[,,t],dvmax=dvmax,nc=nc)    
 trv[,-((t-1)*dvmax+(1:dvmax)),t]<-tr #Assigning values to the related place in the travel array
 }
updtcity<-array(0,dim=c(length(sa[,(nc*dvmax+1),1]),1,nc))  
for(l in 1:nc){
 updtcity[,1,l]<-trv[,nc*dvmax+1,l] }
 for(n in 1:nc){
 sa[,1,n]<-updtcity[,1,n]    
 trv<-trv[,-(nc*dvmax+1),] 
 for(i in 1:(nc-1)){
 for(j in (i+1):nc){
#From i to j sendings
 sa[,1,i]<-sa[,1,i]+sa[,(i-1)*dvmax+2,j]    
 sa[,((i-1)*dvmax+2):((i-1)*dvmax+dvmax),j]<-sa[,((i-1)*dvmax+3):((i-1)*dvmax+dvmax+1),j]  
 sa[,(i-1)*dvmax+dvmax+1,j]<-0
 sa[,(i-1)*dvmax+(1:dvmax)+1,j]<-sa[,(i-1)*dvmax+(1:dvmax)+1,j]+trv[,(j-1)*dvmax+(1:dvmax),i]  
#From j to i sendings
sa[,1,j]<-sa[,1,j]+sa[,(j-1)*dvmax+2,i]
sa[,(j-1)*dvmax+dvmax+1,i]<-0
sa[,(j-1)*dvmax+(1:dvmax)+1,i]<-sa[,(j-1)*dvmax+(1:dvmax)+1,i]+trv[,(i-1)*dvmax+(1:dvmax),j]}}
return(sa)}
generatesend<-function(s,prb,dvmax,nc){ 
#This function returns the matrix includes the number of people travel between each city
   prb<-prb[!is.na(prb)]   
  R<-matrix(0,nrow=dvmax,ncol=length(s))     
  R<-rmultinomv(n=length(s),size=s,prob=prb)
  RN<-t(R)
 return(RN)}
